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ABSTRACT 


Various  procedures  are  given  for  writing  ejtplicit  difference  ap¬ 
proximations  to  the  one-dimensional  Lagrangian  hydrodynamics  equations. 
Computational  comparisons  are  made  among  systems  of  equations  with  tim¬ 
ing  modifications.  These  comparisons  lead  to  experimentally  superior 
differencing  forms.  Stability  analyses  of  these  difference  forms  show 
the  reasons  for  the  superiority  of  one  form  over  another.  Of  greater 
importance,  the  stability  criteria  obtained  show  the  function  of  an  ar¬ 
tificially  introduced  diffusion  term  required  in  the  treatment  here 
given  to  shocks.  The  stability  criterion  in  each  case  involves  the 
familiar  Courant  condition  and  a  term  which  corresponds  to  the  stability 
criterion  of  the  diffusion  equation.  Upper  limits  to  the  magnitude  of 
the  coefficient  of  the  diffusion  term  are  established  as  a  function  of 
Courant  number.  While  lower  limits  are  also  indicated,  they  require 
modification  when  shocks  are  involved. 

Alternate  differencing  schemes  are  considered  in  which  the  pre¬ 
viously-used  total  energy  calculation  is  replaced  by  an  internal  energy 
calculation.  It  is  shown  that  care  must  be  taken  that  the  kinetic  and 
internal  energies  are  expressible  in  terns  of  local  quantities.  That 
is,  in  addition  to  the  equations  being  conservative  in  a  gross  sense, 
they  must  also  be  locally  conservative.  This  is  necessary  in  order  that 
the  energy  condition  of  the  Rankine-Hugoniot  equations  be  satisfied 
when  shocks  arise. 

Finally  discussion  is  given  to  errors  resulting  from  the  replace¬ 
ment  of  shocks  by  a  shock  layer,  that  is,  errors  connected  with  the  ar¬ 
tificially  inserted  diffusion  term.  These  errors  are  manifested  in 
distortions  of  profiles  at  material  discontinuities  through  which  shocks 
have  passed  and  in  rarefactions  associated  with  such  occurrences.  The 
errors  in  turn  effect  stability  in  the  vicinity  of  the  material 
discontinuities. 
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CHAPEER  I 


INTRODUCTION  TO  THE  DIFFERENCING  PROCEDURES 


The  Differential  Equations 


The  difference  approximations  discussed  on  the  following  pages  are 
based  upon  the  Lagrangian  hydrodynamic  equations  written  in  conservative 
form^  with  t,  the  time,  and  m,  the  mass,  as  independent  variables.  The 
mass,  momentvun,  and  energy  equations,  respectively,  are: 


at 


0 


Su  dp 


0 


dE  d(Pu) 
dt 


0 


The  dependent  variables  characterizing  the  fluid  are: 


p  =  density 
u  =  velocity 
P  =  pressure 
E  =  specific  energy 
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The  necessary  additional  equation  is  the  equation  of  state. 

The  coordinate,  x,  of  an  element  of  fluid  in  the  laboratory  freune  is 
a  function  of  m  and  t.  In  terms  of  it  the  velocity  is  defined  by 

ax 

^  =  "St 

and  through  the  mass  equation  it  follows  that 


i  ^  ax 

p  ^ 

Some  discussion  will  be  given  to  the  energy  equation  written  in 
terms  of  the  specific  internal  energy.  This  equation  with  mass  and  time 
as  independent  variables  is 


ai 


0 


where 


Differencing  Procedures 

To  solve  the  initial-condition  boundary-value  problem,  we  replace 
the  differential  equations  by  a  set  of  finite -difference  approximations 
whereby  solution  proceeds  by  a  set  of  purely  algebraic  operations.  To 
accomplish  this  we  represent  the  continuum  of  fluid  by  a  set  of 
"finite  elements"  of  mass  m^,  where  j  =  1,  2,  5  ...  Jq  number  the  cen¬ 
ters  of  mass  of  the  elements  or  "cells."  These  elements  are  shown  in 


the  following  diagram. 


A-J 

_  JL 

j+l 

j-i 


The  Lagrangian  equations  relate  the  field  variables  at  points  mov¬ 
ing  vd.th  the  fluid,  so  that  points  of  information  in  the  mesh  of  cells, 
whether  they  are  centers  of  mass  of  the  elements,  ,3  +  1 

or  cell  boundaries,  j  -  5/2,  j  -  1/2,  j  +  1/2  ••.,  will  include  be¬ 
tween  them  increments  of  mass  which  remain  constant  in  time.  In  the 
following  discussion  we  will  refer  to  the  centers  of  mass  of  the  cells 
simply  as  cell  centers.  It  should  be  noted  that  this  does  not  in  gen¬ 
eral  correspond  to  the  special  cell  center. 

There  is  some  arbitrariness  in  regard  to  the  points  where  the 
field  variables  should  be  defined.  Some  must  be  defined  both  at  cell 
centers  and  at  cell  boundaries.  We  describe  a  procedure  in  which  points 
of  information  and  the  sequence  of  calculations  are  as  follows; 

1 ,  Compute  velocity  changes  at  cell  boundaries  with  the  momentum 
equation, 

2,  Compute  density  changes  at  cell  centers  with  the  mass  equation, 
5.  Compute  energy  changes  at  cell  centers  with  the  energy  equation, 
k.  Compute  cell  center  pressures  with  the  equation  of  state. 

In  the  process  of  differencing  where  the  values  of  the  field  variables 
are  required  and  are  not  given  by  the  above,  we  use  interpolation  for¬ 
mulas  , 

We  tentatively  make  a  direct  correspondence  between  the  differen¬ 
tial  and  the  difference  equations,  that  is,  let 
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(1) 


The  superscripts  n  are  used  to  designate  time  steps*  Thus,  for  example, 
represents  u(nAt),  represents  u[(n+l)At],  etc. 

To  see  more  clearly  what  we  have  done  in  making  the  above  corres¬ 
pondences  of  the  differential  and  difference  equations,  assume  a  know¬ 
ledge  of  the  boxmdary  velocity,  u“  i,  and  the  boundary  pressure, 
at  time  n.  Then  the  other  quantities  in  the  momenttun  equation  in  terms 


of  Taylor  expansions  are 


n+1  n  1^1  j. 

u.  1  =  u.  1  +  Atl'^  j  +  — 5— 
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or  a  first  order  approximation  to  the  derivative  is 


(5) 


Similarly  a  first  order  approximation  to  the  time  derivative  of  the  ve¬ 
locity  is  given  by 


Replacing  the  differential  equation  by  these  approximations  at  the  point 
j-^,n  leads  to 


,n 


n+1 
J  2 


n 

"i-i 
J  2 


p” 

_J - - 


At 


(6) 


which  is  just  the  momentum  equation  (l)  given  above. 

Note  now  that  if  we  add  the  pressures  p'?  and  in  the  above  ex- 

pansions  and  replace  the  mass  derivative  of  the  pressure  as  specified  by 
(5)^  then  we  have 


P^  +  P 

J 


n 

j-1 


=  PP*?  i  + 
0“2 


Thus  a  first  order  approximation  to  the  boundary  pressure  is 


m.  P.  .  +  m.  , 
-  .1  j-' 


,n 


m .  +  m .  , 
J  J-1 


(7) 


This  eqxiation  provides  one  of  the  necessary  intei^olation  formulas  since 
the  cell  boundary  pressures  are  required  in  equation  (4). 

A  similar  treatment  of  the  mass  equation  leads  to  an  Interpolation 


-n 


formula  for  cell  center  velocities,  namely. 


1  +  u.  1 

(8) 

j  2 

Although  the  necessity  for  such  an  interpolation  formula  is  not  indicated 
in  the  above,  we  shall  need  to  evaluate  the  material  pressure  as  a  func¬ 
tion  of  internal  energy,  which  in  turn  is  obtained  by 


In  the  examples  here  discussed  the  polytropic  gas  equation  of  state 
is  used,  i.e., 

P  =  (7  -  1 )pl  (9) 

m 

The  subscript  m  is  used  to  distinguish  the  material  pressure  obtained 
through  the  equation  of  state  from  the  total  pressure  P  which  includes 

p 

a  quantity  q  called  the  "psuedo-viscous  pressure; "  that  is 

P  =  P  +  q  (to) 

m 

The  term  q  is  introduced  into  the  equations  to  spread  shock  fronts 
and  "smooth"  computational  results.  The  spread  of  shock  fronts  to  a 
thickness  on  the  order  of  the  mesh  spacing  avoids  the  necessity  for  spe¬ 
cial  treatment  of  shocks  in  the  difference  approximation.  The  effect  of 
q  must  be  such  that  shock  speed  and  jump  in  the  field  variables  across 
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the  front  are  correct.  This  is  not  in  general  difficult  to  achieve. 

The  Rankine-Hugoniot  jump  condition  equations  are  merely  statements  of 
"local"  conservation  of  mass,  energy,  and  momentum,  that  is,  they  are 
statements  of  the  conservation  laws  across  the  shock.  Hence,  if  the  sys¬ 
tem  of  difference  equations  is  "locally"  conservative,  then  the  Rankine- 
Hugoniot  conditions  are  satisfied  and  the  treatment  will  be  correct  ex¬ 
cept  for  the  spread  of  the  shock  front.  Some  further  discussion  of  this 
will  be  given  later  in  connection  with  con5)utational  experiments. 

Smoothing  of  fluctixations  through  the  use  of  q  has  bearing  upon  the 
stability  of  the  difference  equations.  Considerable  discussion  of  this 
will  be  given  in  Chapter  III.  In  those  cases  where  the  system  of  equa¬ 
tions  is  unconditionally  unstable,  the  addition  of  q  can  make  the  sys¬ 
tem  conditionally  stable.  If  the  system  is  conditionally  stable  without 
q,  the  smoothing  effect  is  still  desirable. 

5 

The  form  of  q  here  considered  is 


n 

j 


C^ 

J 


n 


=  0 


if  u^  1  -  u“  1  >  0 

if  u.  1  -  u.^  <  0 
0“2  d+^  — 


(11) 


C  is  the  material  sound  speed  and  A  is  a  constant,  the  magnitude  of 
which  is  discussed  in  detail  in  Chapter  III.  Note  that  q  is  set  to 
zero  if  the  material  is  in  a  process  of  expanding.  The  desirability  of 
this  will  be  discussed  further  in  Chapter  V  in  connection  with  effects 
upon  rarefactions. 
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Tentative  to  the  discussion  in  Chapter  III,  the  Courant  condition 


C  ^  <  1 
Ax 


will  be  used  as  a  measure  of  stability.  That  is,  we  will  require 


max 


^  <  1 


In  places  where  this  measure  of  stability  is  used  to  control  the  time  in¬ 
crement,  violations  of  the  true  stability  conditions  may  occur,  but  these 
will  be  localized  and  temporary. 
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CHAPTER  II 

COMPUTATIONAL  TESTS  OF  DIFFERENCING  FORMS  INVOLVING 
EVALUATION  OF  TOTAL  ENERGY* 

Prior  to  testing  by  computation  any  particular  form  of  differencing, 
it  is  always  important  to  examine  the  character  of  the  error  terms  aris¬ 
ing  in  the  differencing  procedure.  Usually  it  is  desirable  that  the 
error  terms  vanish  for  infinitesimal  space  and  time  increments.  The 
types  of  differencing  considered  here  all  come  under  this  category  pro¬ 
vided  that  certain  stability  criteria  are  satisfied.  [Note:  For  dis¬ 
cussion  of  a  method  of  accuracy  analysis  the  reader  is  referred  to  Re¬ 
ference  4  by  Harlow  p.  11  ff]. 

In  establishing  what,  in  the  form  of  computational  results,  char¬ 
acterizes  a  good  differencing  scheme,  it  is  helpful  to  have  an  analytic 
solution  of  some  simplified  problem  for  comparison.  If  no  such  solution 
is  known,  an  alternative  basis  for  comparison  is  a  result  obtained  with 


This  is  in  contrast  to  difference  forms  in  which  the  internal  energy 
is  evaluated,  and  the  total  energy  is  treated  as  an  auxiliary  qiian- 
tity.  Such  forms  will  be  discussed  in  Chapter  IV. 
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a  highly  refined  mesh,  that  is,  a  solution  common  to  all  differencing 
forms  under  consideration.  If  the  difference  approximation  represents 
well  some  of  the  more  important  fixnctionals,  then  more  stringent  cri¬ 
teria  may  be  applied.  Practical  considerations  require  that  we  have  as 
much  accuracy  as  possible  in  a  coarse  mesh.  These  considerations  are: 

( 1 )  The  use  of  a  minimum  number  of  space  mesh  points  to  reduce  storage 
requirements  and  thus  also  reduce  calculation  time  per  cycle.  (2)  The 
use  of  a  minimum  number  of  time  cycles  to  reduce  the  over-all  calcula¬ 
tion  time.  Hence  a  loss  of  accuracy  in  a  coarse  mesh  may  be  used  as  a 
device  to  eliminate  some  differencing  forms. 

A  final  computational  criterion  that  may  be  used  is  a  comparison 
in  smoothness  of  results,  that  is,  a  comparison  of  the  rate  at  which 
fluctuations  damp.  If  a  smoothing  mechanism  is  involved,  as  in  the  pre¬ 
sent  discussion,  it  must  be  the  same  for  those  difference  forms  to  be 
compared.  In  using  this  criterion  we  are  basically  looking  for  a  system 
of  difference  equations  which  possesses  error  terms  that  contribute  to 
smoothness  without  affecting  the  accuracy  of  the  over-all  result. 

The  physical  situation  considered  in  the  following  tests  is  that  of 
two  materials;  the  material  on  the  left,  the  energy  source,  is  isother¬ 
mal  and  its  temperature  or  internal  energy  is  specified.  Energy  flows 
from  this  material  into  an  initially  cold,  heavier  material.  A  shock 
wave  moves  into  the  initially  cold  material  and  a  rarefaction  moves  into 
the  isothermal  material. 


First  Test 


We  begin  with  the  system  of  difference  equations  developed  in  Chap¬ 
ter  I  and  concentrate  on  modifications  of  equation  (2).  Equation  (2) 

again  is 


n+1 
X.  1 
0-2 


X .  1  +  u .  1  At 
0-2  0-5 


(12) 


We  compare  it  with 


n+1  n 

X.  1  =  X.  1  + 

0-2  0-2 


n  n+1 

u.  1  +  u.  1 


At 


(15) 


and 


n+1  n  n+1  .. 

X;,  i  =  X  i  +  u  At 
0*2  «3“2 


(1^) 


The  system  of  difference  forms  as  a  whole  with  the  above  modifica¬ 
tions  is  relatively  poor;  our  present  objective  is  only  to  demonstrate 
which  of  (12),  (13)  and  (l4)  is  the  best.  The  results  of  comparisons 
clearly  showed  that  (l4)  was  superior  on  the  basis  of  the  discussed  con¬ 
siderations.  This  was  first  evident  in  calculations  where  the  size  of 
At  was  controlled  by  the  Courant  condition,  that  is,  where  the  time  step 

size  was  adjusted  to  satisfy  C/ax  At  <  1 .  If  for  any  time  step  this 

ins^x 

condition  was  violated.  At  was  reduced  by  half  and  the  calculation  was 
repeated.  On  the  other  hand,  if  we  everywhere  had  C/ax  At  <  5,  then  At 
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was  doubled  on  the  following  time  step.  The  order  of  largest  to  smallest 
time  step  assumed  in  the  three  calculations  was  (l4)^  (15)  (l2). 

Thus  in  (l2)  and  (15)  there  was  a  greater  tendency  toward  fluctuations 
which  in  turn  resulted  in  more  violations  of  the  Courant  condition. 

The  more  coarse  time  mesh  possible  with  (l4)  indicated  that  it  was  the 
more  desirable  form. 

Further  confidence  was  obtained  in  this  conclusion  by  comparing  the 
three  forms  at  the  same  fixed  At.  Typical  velocity  profiles  along  with 
the  theoretical  curve  are  shown  in  Figure  1 .  Details  of  initial  condi¬ 
tions  are  given  in  Table  1 .  The  plot  of  the  rarefaction  to  the  left  of 
cell  1  is  rather  poorly  represented  in  terms  of  Lagrangian  mesh  points 
but  was  of  no  interest  since  it  was  the  same  in  all  cases.  The  relative 
magnitudes  of  the  fluctuations  in  the  region  of  the  initially  cold  mater¬ 
ial  through  which  the  shock  has  passed  (cell  1  and  to  the  right)  clearly 
indicate  that  form  (l4)  has  the  least  fluctuation. 

Note  that  (l2)  is  a  truncation  of  the  Taylor  expansion 


and  (l4)  is  a  truncation  of 
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Note  also  that 


n  n+1 
u  +  u 


implies  a  velocity  at  some  time  between  n  and  n  +  1 ,  If  terms  of  second 
order  in  At  are  omitted, 

u"  +  u“*' 

W  - - a - 

By  a  Taylor  expansion  about  n  +  ^, 


hence  ( 1 5)  is  a  truncation  of 


n+1 
3  2 


=  X 


n 


+  Atl 


3-2 


O(At^) 


We  note  that  the  lowest  order  error  terms  in  (l2)  and  (l4)  are  second 
order,  whereas  in  the  case  of  (l5)  the  lowest  order  error  is  of  third 
order.  From  this  point  of  view  we  conclude  that  form  ( 1 5)  should  have 
given  the  best  results.  The  fact  that  (l4)  gave  the  best  results  can  be 
explained  by  the  second  order  error  having  added  a  contribution  to  the 
smoothing  of  the  results.  A  more  satisfactory  point  of  view,  however,  is 
one  involving  a  relabeling  of  the  equations  such  that  the  time  index  on 
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u  is  everywhere  reduced  by  l/2.  This  requires  a  revision  of  the  initial 
conditions  on  u  and  we  specify  u  ^  =  0.  In  revised  form  the  equations 
(l)»  and  (4)  become 

pH  .  pn 

(15) 


n+p* 

0"2 

„  1 
n*»  0 

J“2 

p*? .  - 

,  ^0.-1.  _._0. 

)  At 

n+1 

X.  1  = 
J-t 

n 

X.  1 

J-2 

+u”*f  At 

pn+1  _ 

{pV-^),,4  - 

(P  U 

0 

+ 

m . 

0 

At 


(i6) 


(17) 


This  revision  in  the  time  index  on  u  now  makes  the  magnitude  of  the  lowest 
order  error  in  (l2),  (l5),  and  (l4)  consistent  with  the  computational  re¬ 
sults*  Note  also  that  the  momentum  equation  (15)  is  now  time-centered  and 
is  correct  to  third  order  in  At* 

The  density  equation  (3)  has  remained  unchanged,  but  the  revised 
time  index  on  u  must  be  taken  into  account  in  the  q  equation  ( 1 1 )  which 
now  becomes 


^  (vt  -  vt) 


q^  =  0 


if  u 


n-^ 

j"2 


u""|^  <  0 

0+4/  - 


(18) 


on  u. 


Compensation  could  be  made  for  the  time  shift  in  initial  conditions 
but  this  is  not  necessary  since  the  effect  is  negligible* 
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Second  Test 


Our  basic  equations  for  further  ejq)erimentation  are  now  those  with 
the  revised  time  index  on  u,  and  we  proceed  by  considering  modifications 
of  equation  (l7)»  We  compare  it  with 


j,n+1 

0 


(P"u  )  i 
_ 


(pV) 


At 


(19) 


and 


/  n  n+|-x  /_n  n+iv 

J  j  m. 


At 


(20) 


n 


where  u  is  now  defined  by 


n+4  n~4 

n  u  +  u  ^ 
u  s - - - 


(21) 


In  all  of  these  cases  we  consider  and  as 


1%  b”  -  ^ 

J  J  2 


n-o 


i"*’  .  e”*' 
0  J 


n+i 


The  same  initial  conditions  of  the  first  test  were  again  used, 
and  the  Courant  condition  was  used  to  control  the  At. 

Figure  2  is  a  plot  of  the  velocity  of  the  interface  between  the 
two  materials.  Note  the  continued  oscillations  in  the  case  of  equation 
(17)  even  at  late  times.  Equations  (l9)  and  (20)  show  rapid  damping  of 
initial  fluctuations,  the  velocity  becoming  the  analytical  value  in  a  short 
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number  of  time  cycles, [For  details  of  the  calculation  refer  to  Table  2], 

On  the  basis  of  the  results  as  demonstrated  by  Figure  2,  eqmtion 
(17)  was  eliminated.  To  contrast  (l9)  and  (20)  more  precisely.  Figure  5 
is  included.  This  is  a  typical,  late-time  density  profile  taken  from  the 
same  calculations.  For  the  present  we  ignore  the  error  in  the  cells  near 
the  material  boundeiry  and  note  that  the  fluctuations  that  are  present 
with  form  (19)  tend  to  average  about  the  profile  of  form  (20).  The  great¬ 
er  smoothness  of  (20)  hence  indicates  that  it  is  the  more  desirable  form 
to  use. 

To  summarize,  consider  the  following  chart  in  which  the  time  in¬ 
dices  on  u  are  used  to  specify  the  types  of  variations  of  differencing 
forms  which  we  have  studied.  The  indices  refer,  respectively,  to  the 
first  and  second  tests. 


^  1 
n-2 

n-2 

n+2# 

n-2 

n-2^ 

n 

n 

n+l, 

n 

„  1 

n+'g- 

n+g- 

n+l. 

n+g 

The  first  test  was  the  upper  row  and  the  second  test  the  right  hand 
column.  Our  conclusions  were  that  the  combination  n  +  §,  n  +  5  was  the 
best  form,  that  being  the  form  where  the  position  equation  and  the 
energy  equation  both  used  the  most  advanced  time  on  the  velocity. 

Four  of  the  nine  variations  in  differencing  listed  above  were  not 
tried  to  this  point.  To  double-check  our  conclusions  the  n,  n  and 
n,  n  +  ^  forms  were  applied  to  the  same  problem.  These  forms  both 
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proved  to  be  poorer  than  both  the  n  +  h  n  and  the  n  +  ^,  n  +  ^  forms. 
This  Indicated  that  there  probably  vas  no  need  to  check  the  remaining 
two  forms^  and  analysis  of  the  next  chapter  confirms  this. 
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CHAPTER  III 


STABIUTY  ANALYSIS 


We  are  interested  in  performing  stability  analyses  of  the  difference 
forms  discussed  in  Chapter  II.  To  include  all  these  forms  we  incorporate 
a  coefficient  n  +  a  on  u  in  the  position  equation,  and  a  coefficient  n  +  p 
on  u  in  the  energy  equation,  where  O!  6ind  p  may  take  on  values  0,  and 

For  simplicity  we  consider  one  material  of  equal  cell  masses.  The  system 
of  equations  \mder  consideration  then  becomes 
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/I 

n  .  „n  f  n-^  n-p 
q.  =  Ap.C.lu  7-“'^  f 
J  \  j-i  j+l:/ 


For  simplicity  in  the  first  analysis  we  are  not  restricting  q  as  in  (ll). 


AssiJine  all  quantities  vary  slightly  from  steady  state  values;  that 


is,  let 


Uj  .  u^d  +  ?  «  1 


Pj  -  Pod  +  €j)  €  «  1 


I.  =  I-(l  +  6  )  6  «  1 

J  0  j 


Then  throu^  the  equation  of  state  we  have  for  first  variations  of  C  and 


c"  =  Vr(r-I)  -  Cgd  +  i  6°) 


(:>„)"  -  (r  -  OpoIo 

u 


The  first  variation  of  q  is  given  by 


n  ,  -  An-4  ^n-4 

'  ‘’o 


The  cell  mass  remains  constant;  hence 


m  =  p6x  =  PQ^^O 


We  also  define  the  Courant  number 


Cq  5t 
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Then  the  first  variations  of  the  momentvun,  mass,  and  energy  ecjuations, 
respectively,  become 


We  now  assume  that  the  solution  of  these  eq,uations  can  be  written 
in  terms  of  a  Fourier  series.  If  this  is  true,  then  each  term  of  the 
series  is  also  a  solution,  and  we  may  examine  a  typical  one  to  note  what 
the  conditions  are  that  would  make  it  a  solution.  Assume  that 


J 


^ikj  n 
e  r 


e  r 


(51) 


Substitution  of  (51 )  into  equations  (28),  (29),  and  (30)  leads  re¬ 
spectively,  to  the  equations 
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(ai  sin|)£  +  jr  [(r  -  1)  +  It  sln^  l]  ^  r  sin  |)  5  -  0 

(r  -  1 )  e  +  ^21  ^  H  sin  ^  =  0  (32) 


u. 


2i(7  -  i)  ^  sin  |j 


C  +  (r  -  1  )  5  =  0 


Then  (31 )  is  a  solution  provided  that 


D  = 


2i  ^  sin  I 
Uq  7  2 


(r  -  1) 


0 


J_ 

/r 


( r  -  1 )  +  4  Tvp.  sin^  ^ 


0  a  .  k 
21  a  r  sin  - 


21(7  -  1)  ^  ^  sin  I 


0  ^  .  k 

2i - sin  - 

Uq  7  2 


0 


(r  -  1) 


=  0 
(33) 


or 

(r 


-  1) 


(r  -  l)  +  4  >s[i  sin^ 


k 


p  p  ih 

4  h  sin  2 


[  w4 


+  (7  -  l)  r' 


0 


(34) 


For  stability  it  is  necessary  that  the  typical  component  solution 
(31 )  not  grow  without  bound.  We  therefore  require  that 


rl  <  1 


(35) 
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We  further  require  that  (35)  satisfied  in  the  extreme  case  of 

p 

sin  k/2  =  1 ;  this  geometrically  implies,  for  example,  that  if  k  =  jt 
the  wave  length  of  the  typical  component  solution  under  consideration 
is  just  twice  the  cell  size.  Under  these  conditions  it  is  possible  for 
this  component  of  each  of  the  field  variables  to  attain  extreme  maximum 
and  minimiim  values  in  alternate  cells,  and  the  largest  possible  gra¬ 
dients  can  exist. 

With  these  considerations  we  wish  to  solve  for  r  the  equation 

(r  -  1)  [(r  -  l)  +  4  Xu]  +  +  (7  "  O  =  0  (36) 

We  shall  confine  our  attention  to  those  values  of  a  and  p  which  lead 
to  a  quadratic  equation  in  r.  These  cases  are  the  ones  which  appear  at 
the  corners  of  the  chart  on  Page  22.  They  are 

a  =  p  =  -  i 

a  =  -  P  =  i 
a  =  i,  p  =  -  i 
a  =  p  =  I 

In  all  of  these  cases  the  equation  reduces  to  a  form 

(r  -  1  )^  +  4^  a(r  -  l)  +  4|i^  =  0  (56a) 

where  a  is  a  function  of  X  and  (i*  The  forms  of  a  for  the  four  cases 
under  consideration  are: 


-28- 


a  =  p  =  -  i 

a  =  A 

a  -  2}  P  -  “  2 

a=?v+^ 

7 

a  -  -  i 

a  =  p  =  i 

a  .2— lJ..|a 

a  =  ^  +  d 

7 

Solved  for  r,  equation  (56a)  becomes 


r 


2d 


a 


(57) 


Now  if  a  <  1 ,  r  is  imaginary  and 


r|  =  Vl  -  4|Ji(a  -  d) 


Thus  the  condition  |r|  <  1  requires  d  <  a. 

If  a  =  1 ,  then  |r|  <  1  requires  0  <  d  <  1 .  Here  the  lower  limit  is 
always  satisfied  while  the  upper  limit  is  included  for  a  <  1 . 

For  a  >  1  we  must  examine  both  the  plus  and  minus  sign  in  the  equa¬ 


tion 


r 


1  -  2da  ±  2d 


1 


(38) 


If  the  plus  sign  satisfies  r  <  1 ,  then  the  minus  sign  will  also.  With 
the  plus  sign,  r  <  1  requires  that 

-  a  +  -  1  <  0 

for  p  >  0.  Note  that  this  is  always  satisfied,  since  a  is  real. 

If  the  minus  sign  in  (38)  satisfies  r  >  -  1,  then  the  plus  sign  will 
also.  Then  r  >  -  1  requires  that 
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This  may  also  be  written  as  a  dual  condition. 


2 

2  a|i  -  h  <1  and  h  <  1 

To  summarize  then,  the  conditions  for  stability  are 


For  a  <  1 ,  <  a 

2 

For  a  >  1  and  ^<l,2ah-U  <1 


(39) 


H- 

a  =  /3  =  -l/2 
a  »X 


X<l  ,  /t<X 

X>l  and  n<\  ,  2X/t-/i?<l 


0=1/2,  ^=-1/2 
o  =  X  +  f 

X+  ■y  >  I  Ond  ^  <  I  ,  2X/i—  5  I 


Plots  of  the  regions  of  stabil¬ 
ity  are  given  for  the  cases  consid¬ 
ered.  Note  that  the  forms  of 
differencing  that  were  experimentally 
superior  are  those  for  which  the  re¬ 
gion  of  stability  is  shifted  to  the 
lowest  values  of  A,  Qualitatively 
the  inference  is  that  the  magnitude 
of  q,  and  hence  the  degree  of  en¬ 
tropy  increase  introduced,  must  com¬ 
pensate  for  certain  of  the  error 
terms  that  are  inherent  in  a  given 
system  of  difference  equations.  If 
the  error  terms  are  of  the  type  that 


cause  an  entropy  decrease,  they  lead 
to  fluctuations  and  hence  require 
more  compensation  in  the  form  of  the 
dissipative  mechanism  q.  [For  fur¬ 
ther  discussion  of  the  relation  be¬ 
tween  the  entropy  of  a  region  and 
fluctuations  the  reader  is  referred 
to  the  Appendix. 3 

It  is  interesting  to  note  that 
the  introduction  of  q  into  the  system 
of  equations  is  equivalent  to  adding 
a  diffusion  term  to  the  momentum 
equation.  That  is, 


a*  -1/2  ,  ;9»l/2  a  =  )9=  1/2 

a  =  X+  fL  0  =  X+^ 

X  +  fj.  <  \  j  fL  <  y\.  X  +  /t5l,(Xj:0) 

X+^^/i.  >  I  8  fL<  I  ,  2X/i  +  ^^^  s  1  X  +  ^>  I  and  I  ,  2X/i.  +  ^^<  I 


St  *  -  S  *  S  ^  35) 

where  ApCbm  is  the  diffusion  coefficient.  The  stability  criterion  for  the 
diffusion  equation  with  this  coefficient  is 

T,  5t  ^  1 

ApCom  — g’  <  ^ 

6m 

which  reduces  to 
2Ap  <  1 


Thus  it  is  seen  that  the  stability  conditions  given  above,  relating  to 
the  upper  limit  on  A,  all  contain  two  terms,  one  connected  with  the 
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hydrodynamic  equations,  the  other  with  the  added  diffusion.  Only  in  the 
case  a  =  p  =  ^  are  the  equations  conditionally  stable  without  any  added 
diffusion,  i.e.,  A  =  0.  It  will  be  noted  that  in  this  case  the  hydrody¬ 
namic  term  alone  is  just  the  Courant  condition. 

In  the  above  analysis  we  did  not  consider  the  dissipative  mechanism 
in  the  sense  specified  in  Chapter  II.  That  is, 

q^  =  0  for  u"^”?  -  u”"i  <  0 

To  extend  the  analysis  to  this  restricted  q,  we  consider  a  typical  cell 
interface  velocity  and  assume  that  the  most  extreme  fluctuations  in 
this  velocity  arise  from  the  existence  of  a  viscous  pressure  to  the  left 
at  one  time  cycle,  and  a  viscous  pressure  to  the  right  on  the  next  time 
cycle.  Let  r^  refer  to  the  times  for  which  the  viscous  pressure  is  from 
the  left,  and  r^  to  the  times  for  which  the  viscous  pressure  is  from  the 
right.  Then  stability  at  time  n  requires  that 


r 


1 


n/2 


<  1 


or 


(40) 


The  equations  of  first  variation  of  6q  for  the  cases  r^  and  r^  become, 
respectively. 


6q 


^  Vo 
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£,q  =  -  q”  =  -  Ap  (t“-i  -  ? 

^  J“2 

The  trial  solution  leads,  respectively,  to 

PqSXo  j-1  Cq  ^ 

.  q  -  -  >*■  r  T  -  '“">5 

pQ^^o  ^  ^0  •/; 

Then  6q  is  the  same  for  r^  and  r^,  and  the  stability  criteria 
|r|  <  1 

again  applies  but  now  with  X  replaced  by  V2»  Thus  the  stability  condi¬ 
tions  for  the  four  forms  of  differencing  discussed  above  become: 


a  =  p  =  -  i 

a=  2 

a  =  i,  p  =  -  1 

X  li 

a  =  77  +  — 

2  7 

a  =  -  i,  P  =  i 

a  =  P  =  1 

A  7-1 

X 

a  =  _  +  1 - |i 

2  7 

a  =  ^  ^ 

with  the  conditions  (59)  on  a.  The  graphical  results  also  hold  with 
the  ordinate  replaced  by  7\/2, 

These  stability  conditions  have  been  verified  computationally,  and 
it  has  been  shown  that  the  stability  of  the  systems  is  very  sensitive 
to  these  conditions. 
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When  shocks  are  involved  the  above  conditions  must  be  satisfied^  but 
also  A  must  be  large  enough  to  prevent  cell  boundaries  from  crossing# 
Thus  far  no  analysis  is  known  in  this  connection  and  we  must  rely  upon 
experimentation#  In  Chapter  V  where  some  additional  adjustments  in  the 
stability  criterion  are  considered,  an  example  indicating  the  character 
of  the  required  lower  limit  on  A  will  be  given# 


CHAPTER  IV 


DIFFERENCING  PROCEDURES  INVOLVING  THE  INTERNAL  ENERGY  CALCULATION 


In  practical  applications  of  the  Lagrangian  difference  equations, 
there  has  been  some  interest  in  the  use  of  the  energy  equation  in  the 
form  which  gives  the  internal  energy  directly.  The  reason  is  the  economy 
of  calculations  that  results.  Rather  than  evaluating  both  total  and 
kinetic  energies  to  obtain  the  internal  energy,  it  is  obtained  directly. 
An  additional  advantage  is  that  in  some  formulations  the  interpolations 
need  not  be  carried  out  in  the  calculation  procedure.  Thus  a  worthwhile 
reduction  in  the  number  of  mathematical  manipulations  is  realized. 

We  proceed  with  a  test  of  the  system  of  equations  (l5)>  and 

(5)>  with  the  energy  eq^uation  (l7)  replaced  by  an  internal  energy  calcu¬ 
lation.  The  differential  equation  and  a  first  possible  difference  form 
are 


hi 

'St  =  " 


n 


m . 
0 


Solved  for  the  internal  energy  at  the  advanced  time,  this  is 


in+1 

J 


i^  +  p‘: 
j  j 


n  n 

u.  1  -  u.  1 

^  ^ _ m 


m , 


At 


(41) 


-35- 


where  again 


n 


u 


n- 

+  u 


1. 

2 


2 


We  now  compare  ( if  1 )  with 


1- 
n“f  o 


jU+l 


=  l"  + 

J  J 


r  o  n+ o’ 

i  -  ^..1 

•1~2  “ 


m. 


^  At 


(42) 


From  the  conclusions  reached  in  the  preceding  chapters  regarding 
the  time  index  on  we  exclude  the  case  n  - 

For  computational  comparisons  of  these  two  forms  of  the  energy  equa¬ 
tion^  the  same  initial  conditions  were  again  used*  [For  details  of  the 
calculation  refer  to  Table  2].  Typical  profiles  of  the  internal  energy- 
in  the  initially- cold  material  are  given  in  Figure  4.  The  analytic  re¬ 
sult  is  included  for  conqoarison.  Form  (42)  gave  results  which  are  ser¬ 
iously  in  error,  whereas  form  (4l )  matched  the  analytic  solution  except 
for  fluctuations.  Examination  of  the  other  field  variables  likewise 


It  will  be  noted  upon  examination  of  the  character  of  equations 
(41)  and  (42)  that  if  a  cell  j  is  initially  cold,  it  will  remain  so 
indefinitely  with  the  form  of  q  of  equation  (l8).  To  perform  tests  of 
these  equations  it  was  necessary  to  modify  q  to  permit  starting  condi¬ 
tions.  This  was  done  here  by  replacing  in  q  by 


_J _ 


n 

^.1  “  ^0 


where  u„  is  the  velocity  ahead  of  the  shock  layer.  This  modification 
does  not  affect  the  qualitative  feature  of  the  results  in  terms  of  the 
discussion  here  given.  The  predominant  effects  of  this  modified  q  are 
that  starting  conditions  are  provided  or  augmented,  and  additional 
dan5)ing  of  shock-produced  fluctuations  is  realized. 
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showed  correct  results  with  (4l )  and  incorrect  results  with  (42).  In 
particular  it  is  evident  in  late  time  plots  that  the  shock  speed  in  the 
case  of  (42)  is  lower  than  the  theoretical  value. 

The  reason  for  the  poor  behavior  of  (42)  may  best  be  understood  by 
an  examination  of  the  corresponding  total  and  kinetic  energy  expressions. 
It  can  be  shown  that  fonn  (4l),  which  gave  correct  results,  is  derivable 
from  the  energy  equation  (l9)  with  the  momentiim  equation  (15)  and  inter¬ 
polation  formula  (y)*  It  is  required,  however,  that  the  cell  specific 
kinetic  energy  be  defined  by 


(45) 


rather  than  through  the  interpolation  formula  (8);  that  is,  the  cell 
center  velocity  is  defined  through  an  average  of  the  cell  boundary  ki¬ 
netic  energies  rather  than  an  average  of  the  boundary  velocities.  Of 
importance,  however,  is  that  the  kinetic  energy  of  the  cell  depends 
only  upon  the  cell  boundary  velocities  at  the  particular  time  in  ques¬ 
tion.  By  contrast  an  analogous  derivation  of  equation  (42)  is  possible 
using  equations  (20),  (15)^  and  (7))  but  in  this  case  the  required  de¬ 
fining  formula  for  the  kinetic  energy  is 


n*=1 


(Wt) 
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where  the  index  1 /2  on  the  first  term  on  the  right  is  a  time  index.  Note 
that  this  definition  of  the  cell  kinetic  energy  is  in  terms  of  the  cell 
boundary  velocities  but  involves  all  prior  times  to  the  time  of  interest. 

To  analyze  the  significance  of  the  difference  in  the  above  required 
definitions  of  the  kinetic  energy,  we  diverge  briefly  to  consider  the 
meaning  of  conservation  in  a  finite  difference  scheme.  Since  we  are  con¬ 
sidering  a  net  of  points  of  information,  the  idea  of  conservation  may  be 
resolved  to  the  requirement  that  the  value  of  flux  of  momentum  euid 
energy  remain  the  same  whether  observed  in  connection  with  the  cell  to 
the  left  or  the  cell  to  the  right.  Likewise  the  value  of  a  given  con¬ 
servative  variable  specified  at  half-integer  times  must  remain  the  same 
whether  observed  in  connection  with  time  earlier  or  later  than  the  spe¬ 
cified  time.  That  is,  conservation  is  equivalent  to  the  require¬ 
ment  that  single-valuedness  of  the  variables  exist  at  all  points  of  the 
mesh.  In  our  present  discussion  the  temporal  criterion  is  satisfied  in 
regard  to  the  kinetic  energy,  since  it  is  uniquely  defined  at  the  net 
point  j,n  +  f  by  equation  (44).  The  problem,  however,  is  that  when 
shocks  are  involved,  the  Rankine-Hugoniot  equations  must  be  satisfied 
for  the  jump  in  the  field  variables  at  the  shock.  The  Rankine-Hugoniot 
equations  are  statements  of  conservation  of  mass,  momentum,  and  energy 
across  the  shock  front.  Since  the  shock  is  a  local  variation,  these 
equations  require  local  conservation.  Local  conservation  of  energy  re¬ 
quires  definitions  of  kinetic  and  internal  energies  in  terms  of  local 
quantities.  Thus  in  equation  (44)  the  summed  terms  should  not  be 
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present.  We  note  that  the  difference  in  the  kinetic  energies  at  two  ad¬ 


jacent  times  is  given  by 


Hence  in  a  local  sense  the  kinetic  energy  as  observed  from  time  n  -  1  at 
time  n  -  5  is  given  by 


and  as  observed  from  time  n  is  given  by 


2  2 


Thus  it  is  seen  that  locally  the  kinetic  energy  is  not  single-valued 
(that  local  conservation  does  not  exist),  and  hence  the  Rankine-Hugoniot 
equations  are  not  satisfied. 

From  this  discussion  we  see  that  it  is  important  that  the  form  of 
the  finite  difference  expressions  be  founded  upon  local  definitions  of 
those  quantities  which  must  be  conserved.  This  can  be  achieved  in  the 
present  discussion  by  beginning  with  a  system  of  difference  equations 
which  involves  a  total  energy  calculation.  In  such  a  case  an  explicit 
definition  of  the  kinetic  energy  must  be  made,  and  hence  the  require¬ 
ment  of  local  definition  may  be  imposed.  Such  a  system  of  equations  may 
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then  be  used  in  deriving  short-cut  procediures. 


Although  no  problems  of  the  type  discussed  above  arose  in  the  pre¬ 
sent  formulation  of  the  momentum  equation,  the  same  considerations  apply 
to  it  in  general. 

We  now  wish  to  examine  some  examples  of  internal  energy  calculations 
based  upon  the  total  energy  equations  given  in  the  foregoing  discussion. 
There  are  a  nvunber  of  possibilities,  but  we  will  consider  only  four  of 
them.  These  four  differ  in  the  time  index  on  u  and  in  the  kinetic  energy 
interpolation  used.  We  may  tabulate  these  cases  with  n  or  n  +  ^  specify¬ 
ing  the  time  index  on  u  in  the  total  energy  equation,  and  (8)  or  (43) 
specifying  the  kinetic  energy  interpolation  formula.  Thus  we  designate 


(a) 

(8) 

(b) 

(43) 

(c) 

n  +  1, 

(8) 

(d) 

n  +  i. 

(43) 

In  all  cases  the  boundary -pressure  interpolation  formula  is  equation  (?), 
and  the  time  index  n  +  ^  on  u  is  used  in  the  position  computation. 

The  derived  internal  energy  equations  are 


(a) 

^n+1 

H 

il 

At 

/„n  nv  /_,n  n>. 

1  (P  u  )^_|  -(Pa) 

^n+1 

At 

r/  n  ,  n  A 

+  — — 
m . 

1 

'2 
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(c) 

_n+1 

J 

H 

II 

(pV^2)  ^  . 

L  0-2 

/„n  n+^v  n 

- 

r ,  .2 

^  N. 

(d) 

_n+1 

J 

-3 

m . 

J 

n/  n-|^ 

Plu  ,  -  u  f 
dV  J-2  J+2-^ 

) 

fe-i)  *  ' 

Gm)  ] 

At^ 

The  index  n  on  u  wherever  it  appears  again  means 


^  n+4  n-4 

n  u  +  u 
u  =  - p; - 


In  (d)  the  acceleration  a^  is  given  by 

J~2 


n  ^J-1  ■ 


We  have  introduced  only  one  new  difference  form,  (d).  (a)  and  (c) 

were  discussed  in  Chapter  II  where  they  appeared  as  equations  (l9)  and 
(20),  respectively,  (b)  is  equation  (4l )  repeated  here  for  comparison. 

Note  that  the  new  form  (d)  without  the  last  term  is  Just  equation 
(42).  The  second  order  term  in  At  added  to  (42)  resulted  from  our  de¬ 
finition  of  the  kinetic  energy  in  terms  of  local  quantities  (equation 

43). 

Of  passing  interest  we  note  that  (a)  and  (c)  are  difference  approxi¬ 
mations  to 

bl  _  r^(Pu)  dP 

^  ^  ^ 


They  differ  only  in  the  time  index  on  u  in  the  work  flux  at  the  cell 
boundaries.  In  Chapter  II  we  concluded  that  (c)  was  a  superior  form  to 
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use  in  that  it  led  to  smoother  profiles  for  a  given  q  and  At  (refer  to 
Figure  3).  Likewise,  while  the  stability  conditions  for  (a)  were  not 
given  in  Chapter  III,  the  trend  indicated  that  the  most  advanced  time 
required  the  least  compensation  of  errors  in  terms  of  the  magnitude  of  q. 
The  stability  conditions  for  (c)  and  (d)  are  the  same  and  are  the 
a  =  p  =_l/2  case  discussed  in  Chapter  III.  (a)  and  (b)  correspond  to 
a  =  1 /2,  p  =  0  in  the  notation  there  used.  Differences  in  the  form  of 
the  kinetic  energy  interpolation  foimula  vanish  when  the  equations  are 
linearized,  thus  leading  to  identical  stability  conditions  through  the 
type  of  analysis  of  Chapter  III. 

To  further  check  the  assiimption  that  the  more  advanced  time  on  u 
gives  superior  results,  forms  (b)  and  (d)  were  compared  using  the  same 
initial  conditions  of  previous  tests.  [For  details  of  this  test  refer 
to  Table  2].  For  this  comparison  typical  density  profiles  are  given  in 
Figure  5  along  with  the  analytic  result.  Again  we  ignore  the  error  in 
the  cells  near  the  energy  source  material  and  make  the  comparison  on  re¬ 
lative  smoothness  of  the  profiles.  We  note  that  (d)  tends  to  take  on 
mean  values  of  (b)  and  is  therefore  considered  to  be  superior. 

To  take  further  advantage  of  the  advanced  time  on  u,  it  might  be 

n+— 

suggested  that  u^  in  form  (c)  be  replaced  by  u"  2.  This,  however,  leads 
back  to  the  problem  of  the  kinetic  energy  being  defined  in  terms  of  non¬ 
local  quantities. 

A  comparison  of  (c)  and  (d)  is  of  some  interest.  Although  these 
forms  are  the  same  in  terms  of  first  variations,  they  differ  in  second 


order.  Form  (c)  differs  from  (d)  in  that  the  kinetic  energy  is  always 
less  by  a  factor 


This  reduction  in  kinetic  energy  goes  into  internal  energy  and  hence  re¬ 
sults  in  an  entropy  increase.  This  term  has  most  significance  where  the 
velocity  gradients  are  Isorge  and  thus  acts  as  an  auxiliary  viscosity  in 
shock  regions;  as  would  be  expected,  it  helps  to  prevent  cell  boundaries 
from  crossing  as  a  result  of  shocks. 

It  will  be  noted  that  in  an  adiabatic  rarefaction  the  term  will 
cause  an  entropy  increase  since  it  will  remain  positive.  Although  this 
effect  may  be  small,  it  is  still  undesirable.  It  is  therefore  recommended 
that  if  the  characteristics  of  form  (c)  are  otherwise  desirable,  the 
above  term  should  be  taken  into  account  only  in  cells  for  which 

u.  1  i>0,  as  has  been  specified  for  q.  This  is  equivalent  to 

J“2  0+2 

using  form  (d)  normally,  but  replacing  it  by  form  (c)  in  compressive  re¬ 
gions.  The  merits  of  form  (c)  may  perhaps  more  easily  be  realized  by 

3 

modifying  q  as  suggested  by  Landshoff,  that  is,  by  adding  a  term  to  q 
involving  (Au]P .  Doing  this  will  provide  the  same  effect  in  isothermal 
regions  as  well. 
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CHAPTEB  V 

DISCUSSION  OF  ERRORS 

Whereas  the  difference  methods  presented  in  the  previous  chapters 
give  satisfactory  results  in  most  problem  situations,  there  are  errors 
evident  in  all  cases  which  stem  from  the  artificially  introduced  diffu¬ 
sion  term,  i.e.,  from  the  use  of  the  mechanism  q.  We  will  first  note  the 
form  these  errors  take  in  some  sample  results. 

In  Figure  6  a  density  profile  of  differencing  form  (c)  is  given. 
Reference  to  Table  3  shows  that  this  test  was  made  with  equal  mass  zones. 
That  is,  the  light  isothermal  material  had  large  cells,  while  the  ini¬ 
tially  cold,  heavy  material  had  relatively  small  cells.  The  tests  were 
made  in  this  way  to  minimize  the  truncation  error  in  the  mass  derivatives 
at  the  material  interface.  In  contrast,  a  corresponding  test  (c-4)  was 
made  in  which  cell  sizes  were  more  of  the  same  order  [refer  to  Table  4]. 
The  density  profile  of  this  test  is  also  given  in  Figure  6. 

The  phenomenon  of  interest  is  the  low  density  relative  to  the  theo¬ 
retical  value  that  occurs  in  the  first  few  cells  of  the  shocked  material. 
In  the  case  of  the  cells  of  equal  mass,  the  effect  is  spread  over  several 


cells,  while  in  the  case  of  cells  more  equal  in  physical  size,  the  effect 
is  more  severe  but  confined  to  the  first  two  cells. 

In  Figure  7  "the  internal  energy  profiles  of  the  shocked  material  are 
given  for  the  same  tests.  Note  that  the  effect  under  consideration  is 
here  manifested  in  an  excess  of  internal  energy  in  the  first  cells. 

In  Figure  8  velocity  profiles  in  the  rarefaction  region  of  the  iso¬ 
thermal  material  are  presented  along  with  the  theoretical  curve.  [Refer 
to  Table  5  for  details  of  this  test].  The  solid-line  curve  corresponds 
to  tests  of  difference  form  (c)  with  q  restricted  to  compressed  cells, 
while  the  dashed  curve  is  a  result  in  which  q  was  unrestricted,  i.e.,  q 
took  on  both  positive  and  negative  values.  We  first  note  that  both  ex¬ 
perimental  curves  lag  behind  the  theoretical  curve.  This,  however,  is 
of  no  concern  since  the  lag  is  related  to  starting  conditions  and  does 
not  grow  in  time.  We  note,  however,  that  the  unrestricted  q  causes  an 
unrealistic  spread  in  the  rarefaction  front.  This  incidentally  is  the 
reason  that  the  restricted  q  is  generally  used.  Further,  it  is  noted 
that  while  the  restricted  q  curve  has  the  correct  velocity  gradient  in 
the  rarefaction  fan,  there  is  a  rise  in  the  velocity  immediately  behind 
the  rarefaction.  This  is  not  a  consequence  of  q  in  the  vicinity  of  the 
rarefaction  fan;  in  fact,  q  is  instrumental  in  reducing  a  larger  rise  in 
velocity  in  this  region  that  would  be  present  without  the  viscosity. 

We  must  therefore  relate  this  effect  with  errors  that  occurred  at  the 
origin  of  the  disturbance. 

Finally,  Figure  9  shows  a  density  profile  of  the  differencing  types 
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(c)  and  (d)  of  Chapter  IV  applied  to  a  shock  moving  into  a  material  dis¬ 
continuity  where  the  material  on  the  left  is  more  dense  by  a  factor  of 
25  [refer  to  Table  6  for  details  of  this  test].  When  the  shock  strikes 
the  material  discontinuity,  a  shock  moves  into  the  light  material  on  the 
right  and  a  rarefaction  moves  into  the  more  dense  material  on  the  left. 

The  errors  that  are  here  evident,  relative  to  the  theoretical  curve,  may 
be  correlated  to  the  effects  pointed  out  above.  Again  the  density  in 
the  first  cells  of  the  shocked  material  is  low.  The  excessive  velocity 
behind  the  rarefaction  fan  is  manifested  as  a  region  of  lower  density 
than  the  theoretical  value.  The  high  peak  in  the  density  at  the  mater¬ 
ial  interface  in  the  rarefied  material  is  merely  a  consequence  of  the 
oversized  cells  to  the  right. 

From  the  above  examples  and  discussion  it  is  apparent  that  the  major 

errors  may  be  correlated  with  effects  occurring  at  material  discontinuities. 

The  low  density  and  high  internal  energy  occurring  in  the  first  cells  of 

5 

the  shocked  material  is  discussed  briefly  by  Landshoff.  This  effect  is 
attributed  to  an  overproduction  of  entropy  in  the  affected  cells,  result¬ 
ing  in  excessive  heating  and  a  consequent  ejqiansion  of  these  cells.  It  is 
therefore  related  to  the  character  of  the  dissipative  mechanism  q.  Since 
the  peculiarities  connected  with  rarefactions  were  traced  back  to  where 
the  rarefaction  originated,  then  these  errors  likewise  are  assumed  to  be 
a  consequence  of  the  action  of  q  at  material  discontinuities. 

Although  an  artificial  viscosity  is  essential  to  the  type  of  treat¬ 
ment  here  given  to  shocks,  its  secondary  effects  can  become  of  some 
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concern.  Basically,  q  was  mauiufactured  to  effect  a  replacement  of  a 
steady  state  shock  by  a  shock  layer  of  thickness  on  the  order  of  a  cell 
size.  The  process  of  development  or  modification  of  this  layer  is  not 
provided  for  in  the  character  of  q.  The  errors  evident  in  Figures  6 
and  7  are  correlated  with  the  development  of  the  shock  layer  from  ini¬ 
tial  conditions^  while  in  Figure  9  the  shock  layer  required  modification 
at  the  material  interface. 

Experiments  in  which  the  viscosity  coefficient  A  was  changed  show 
that  a  larger  A  increases  the  entropy  production  at  material  discontinui¬ 
ties  and  hence  increases  the  errors.  Modification  of  A  alone,  however, 
is  not  sufficiently  effective  to  eliminate  the  errors.  A  balance  must 
be  achieved  in  the  size  of  A,  in  which  the  further  considerations  of 
fluctuations  and  crossing  of  boundaries  must  be  considered. 

In  Chapter  III  it  was  pointed  out  that  the  stability  conditions  of 
the  equations  required  modification  when  shocks  were  involved;  A  must 
be  sufficiently  large  to  prevent  boundaries  from  crossing.  An  additional 
modification  is  required  because  of  the  secondary  effects  of  q  discussed 
above. 

In  Figure  10  experimental  points  are  plotted  for  difference  form 
(d)  with  the  initial  conditions  of  Table  7*  These  points  are  superim¬ 
posed  upon  the  previously  discussed  plot  of  the  region  of  stability  for 
difference  forms  a  =  P  =  1  /2,  in  which  q  is  the  restricted  form.  The 
results  must  be  teiken  only  in  a  qualitative  sense,  since  boundary  cross¬ 
ing  is  strongly  dependent  upon  initial  conditions.  The  points  of 
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instability  resulting  from  shocks  are  classified  into  two  groups.  Boun¬ 
dary  crossings  occurred  after  only  a  few  time  steps  and  involved  the 
interface  between  the  isothermal  material  and  the  initially  cold  material. 
Instabilities  arising  from  the  overproduction  of  entropy  occurred  at  re¬ 
latively  late  times  and  were  effective  in  the  region  directly  following 
the  cells  in  which  the  low  density  errors  occurred.  This  latter  type  of 
instability  depends  upon  the  relative  cell  sizes  of  the  two  materials  of 
the  interface.  If  the  cells  of  the  isothermal  material  were  larger,  this 
type  of  instability  would  be  less  likely  to  arise  because  the  low  density 
errors  would  be  less  severe  [refer  to  Figures  6  and  73 •  Note  that  for 
the  type  of  instability  considered,  improvement  with  reduced  ^  is  small. 
The  present  point  of  view  is  that  the  secondary  effects  of  q,  are  such 
as  to  modify  the  condition 


which  certainly  must  be  satisfied  for  compressed  cells. 

In  Chapter  III  it  was  noted  that  a  shift  in  the  region  of  stability 
to  lower  values  of  7\  was  a  qmlitative  indication  that  the  difference 
form  was  superior.  The  question  might  be  raised,  why  then  use  a  re¬ 
stricted  q  which  shifts  the  region  of  stability  to  A's  of  twice  the 
value  of  the  unrestricted  q?  If  a  technique  could  be  devised  to  prevent 
the  appearance  of  q  in  rarefactions,  but  still  retain  it  in  unrestricted 
form  in  all  other  regions,  then  an  improvement  could  be  realized. 
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TABLE  1 


INITIAL  CONDITIONS  FOR  CALCULATIONS  SHOWN  IN  FIG.  1 

ist  fixed  at  O.0625 


Material  No. _ 1  ( isothermal) _ 2 


7 

5/5 

5/3 

1 

1 

0.1 

1 

I. 

J 

0.246 

0 

X 

1/7 

1/7 
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Fig.  1 :  First  test  of  Chapter  II.  Cell  velocity  profile  at  time  t  =  15* 
Cells  0,  -1,  -2,  etc.,  axe  larger  than  cells  to  right  by  a  fac¬ 
tor  of  1 0.  [Refer  to  Table  1 ] . 
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TABLE  2 


INITIAL  CONDITIONS  FOR  CALCULATIONS  SHOWN  IN  FIGS.  2,  3,  4,  and  5 


At  adjusted  to  satisfy 


At  <  1 


Material  No, 


1  (isothermal) 


0.2i<-6 


Refer  to  footnote  Page  56  in  connection  with  Figs.  4  and  5. 
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0  5  10  15  18 


TIME 

Second  test  of  Chapter  II.  Material  interface  velocity  as  a 
fimction  of  time.  [Refer  to  Table  2l. 


CELL  NUMBERS 


Chapter  IV  test.  Internal  energy  profile  of  material  2  at  time 
t  =  15«  [Refer  to  Table  2]. 


CELL  NUMBERS 


Fig.  5:  Chapter  IV  test.  Density  profile  of  material  2  at  time  t 
[Refer  to  Table  2]. 


TABLE  5 


INITIAL  CONDITIONS  FOR  CALCULATION  (c)  SHOWN  IN  FIGS.  6  AND  7 


At  adjusted  to  satisfy  46  <  1 


TABLE  4 


INITIAL  CONDITIONS  FOR  CALCULATION  (c-4)  SHOWN  IN  FIGS.  6  AND  7 


Material  No.  1  (isothermal)  2 


7 

5/5 

5/3 

“.3 

0.i4 

1 

"o 

0.1 

1 

0.246 

0 

A 

1/7 

1 .7 
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CELL  NUMBERS 


Chapter  V  error  illustration.  Density  profiles  of  material  2 
taken  for  different  initial  cell  sizes  in  material  1  ft  =  1 5) 
[Refer  to  Tables  3  and  4]. 


5 


10 


15 


CELL  NUMBERS 

7:  Chapter  V  error  illustration.  Internal  energy  profiles  of 
material  2  taken  for  different  initial  cell  sizes  in 
material  1  (t  =  15)*  [Refer  to  Tables  3  and  k]. 
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TABLE  5 

INITIAL  CONDITIONS  FOR  CALCULATIONS  SHOWN  IN  FIG.  8 


At  adjusted  to  satisfy  *v/6  At  <  1 

^  '^max 

Material  No. _ 1  (isothermal) _ 2 


y  5/3  5/5 

0.25  1 

I  0.246  0 

J 

A  1/7  i/r 
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VELOCITY 


0.12 


POSITION 


Fig.  8:  Chapter  V  error  illustration.  Velocity  profiles  in  material  1 
comparing  two  forms  of  q  in  the  rarefaction  (t  =  15)»  [Refer  to 
Table  5 ] • 
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TABLE  6 


INITIAL  CONDITIONS  FOR 

CALCULATIONS  SHOWN  IN  FIG.  9 

^  adjusted  to  satisfy  46 

'^■4i60C 

/ist  <  1 

Material  No. 

1  (isothermal) 

2.  _ 

5 

7 

5/5 

5/5 

5/5 

m . 

0 

0.25 

1 

0.2 

"j 

0.1 

1 

.04 

^j 

0.246 

0 

0 

A 

i/r 

1/7 

1/7 

*The  form  of  q  described  in  the  footnote  of  page  36  was 
used  here.  It  should  he  noted,  however,  that  neither 
form  (c)  nor  (d)  reqviire  starting  conditions. 
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CELL  NUMBERS 


Chapter  V  error  illustration.  Density  profiles  at  time 
t  =  t^  +  1^.9 f  where  t^  is  time  shock  passed  the  second 
material  interface.  [Refer  to  Table  61. 


TABLE  7 

INITIAL  CONDITIONS  FOR  CALCUIATIONS  USED  IN  CONNECTION  WITH  FIG.  10 
A  and  n  (fixed  At)  as  indicated  on  Fig.  10 


Material  No. 


1  (isothermal) 


m . 


I. 

J 


5/3 

0.25 

0.1 

0.246 


5/3 
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Fig.  10:  Chapter  V.  Example  of  experimental  points  of  stability  and 
nonstability  when  shocks  are  involved.  [Refer  to  Table  7]* 


APPENDIX 


ENTROPY  OF  A  REGION  OF  FLUID 

Consider  a  one -dimensional  region  of  a  fluid  which  is  subject  to 
variations,  the  variations  in  turn  being  subject  to  the  constraints  of 
fixed  mass,  fixed  momentum,  and  fixed  energy.  That  is, 
b 

m  =  /a,  p  dx  =  constant 
b  — 

M  =  /^  p  udx  =  /a  M(p  ,u)dx  =  constant  (A-1 ) 

^  12  b  — . 

E  =  /a  p[l(p,P)  +  2  u  ]  dx  =  /g^  E(p,P,u)dx  =  constant 

where  M  and  E  are  momentum  and  energy  densities,  respectively. 

Now  €ilso  let  S(p,P)  be  the  entropy  density  of  the  fluid  region 
such  that  the  total  entropy  is  given  by 

S  =  /a  S  dx  (A-2) 

We  ajre  interested  in  examining  the  character  of  extremal  values  of 
S  in  terms  of  variations  of  p,  P,  and  u  which  are  arbitrary  except  for 
the  specified  conservative  constraints.  The  condition  for  an  extremal 
value  of  S  is 
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‘fa.  (If  '"P  ^  ° 


(A-3) 


T?ie  constraints  in  terms  of  variations  of  p,  u,  and  P  are 


6p  dx  =  0 


Ja  W  SI  6“  )  to  =  0 


(A-4) 


^  ^  dE  ,  bE 


/;( 


^  6p  .  ^  6u  .  5P  )  dx  =  0 


Multiplying  these  constraints  successively  by  the  Lagrangian  multipliers 
f,  g,  and  h  and  adding  to  the  equation  dS  =  0,  we  have 


(A-5) 


Since  all  the  specified  constraints  are  now  included,  the  variations 
6p,  6P,  and  6u  are  arbitrary  (i.e.  independent);  hence  the  coefficients 
of  these  variations  must  vanish  individually.  We  thus  have 


Ss  ^  v, 

dS  ,  S  . 


du 


g  -  +  0 


dE 


(A-6) 
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Through  specification  of  material  equations  of  state  and  the  above 
definitions,  S,  M,  sind  E  become  known  functions  of  p,  P,  and  u.  Hence 
the  equations  (A-6)  reduce  to  a  set  of  algebraic  equations  in  the  three 
unknowns  p,  P,  and  u  with  solutions 

p  =  constant 
P  =  constant 
u  =  constant 


With  the  specified  constraints  governing  the  behavior  of  the  region  of 
fluid  under  consideration,  we  thus  conclude  that  the  entropy  will  have 
an  extremal  value  characterized  by  a  flat  profile  in  all  variables. 

Since  we  have  established  what  characterizes  this  extremal  value  of 
S,  the  examination  of  any  other  condition  will  tell  us  if  the  extremal 
is  a  maximum  or  a  minimum. 

Consider  a  small  space  variation  in  internal  energy  only,  that  is, 

I  =  Iq  +  €(x) 

P  -  Pq  (A-7) 


Then  the  fixed  energy  constraint  becomes 


f 


T  \  12 

Iq  2  “o 


dx  =  p, 


T  1  2 

^0  *  2  “o 


]/: 


e(x)  dx  =  0 


or 


(A-8) 
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The  remaining  two  constraints  are  automatically  satisfied  since 
(A-7)  specifies  no  variation  in  p  and  u.  Also  for  this  special  case 
S  =  S(l)  only,  and 


With  equation  (A-8)  taken  into  account,  the  change  in  entropy  of 
the  region  becomes 


dS 


S  -  Sq)  dx  «  ^ 


e(x)^  dx 


Note  now  that  the  established  extremvun  of  S  is  a  maximum  if 

is  negative  definite.  With  the  density  constant  the  second 
law  of  thermodynamics  states  that 


is  a  fundamental  property  of  actual  media,  and  therefore 


<  0 
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Thus  we  have  shown  that  a  flat  profile  of  the  variables  character¬ 
izing  the  fluid  implies  an  entropy  maximum  if  the  conservation  laws  hold. 

Conversely,  if  we  consider  a  true  solution  to  the  hydrodynamic  equa¬ 
tions  characterized  by  some  value  of  entropy  and  compare  with  it  a  finite 
difference  approximation  solution,  then  the  degree  of  fluctuation  in  the 
approximation  may  be  correlated  with  error  terms  contributing  to  a  de¬ 
crease  in  entropy.  This  is  assuming  that  the  difference  approximation 
is  conservative  in  the  above  specified  quantities. 
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